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Abstract 

We present asymptotically faster approximation algorithms for the generalized flow prob- 
lems in which multipliers on edges are at most 1. For this lossy version of the maximum 
generalized flow problem, we obtain an additive e approximation of the maximum flow in time 
0(m 3 / 2 log 2 (C//e)), where m is the number of edges in the graph, all capacities are integers in 
the range {1, . . . , U}, and all loss multipliers are ratios of integers in this range. For minimum 
cost lossy generalized flow with costs in the range {1, . . . , U}, we obtain a flow that has value 
within an additive e of the maximum value and cost at most the optimal cost. In many param- 
eter ranges, these algorithms improve over the previously fastest algorithms for the generalized 
maximum flow problem by a factor of m 1 / 2 and for the minimum cost generalized flow problem 
by a factor of approximately m 1 / 2 /e 2 . 

The algorithms work by accelerating traditional interior point algorithms by quickly solv- 
ing the linear equations that arise in each step. The contributions of this paper are twofold. 
First, we analyze the performance of interior point algorithms with approximate linear system 
solvers. This analysis alone provides an algorithm for the standard minimum cost flow problem 
that runs in time O (m 3 / 2 log 2 U) — an improvement of approximately O (n/m 1 ^ 2 ) over previous 
algorithms. 

Second, we examine the linear equations that arise when using an interior point algorithm 
to solve generalized flow problems. We observe that these belong to the family of symmetric 
M-matrices, and we then develop O (m)-time algorithms for solving linear systems in these 
matrices. These algorithms reduce the problem of solving a linear system in a symmetric M- 
matrix to that of solving O (logn) linear systems in symmetric diagonally-dominant matrices, 
which we can do in time O (m) using the algorithm of Spielman and Teng. 

All of our algorithms operate on numbers of bit length at most O (lognU /e). 



*This material is based upon work supported by the National Science Foundation under Grant Nos. CCF-0707522 
and CCF-0634957. Any opinions, findings, and conclusions or recommendations expressed in this material are those 
of the author(s) and do not necessarily reflect the views of the National Science Foundation. 



1 Introduction 



Interior-point algorithms are one of the most popular ways of solving linear programs. These 
algorithms are iterative, and their complexity is dominated by the cost of solving a system of linear 
equations at each iteration. Typical complexity analyses of interior point algorithms apply worst- 
case bounds on the running time of linear equations solvers. However, in most applications the 
linear equations that arise are quite special and may be solved by faster algorithms. Each family 
of optimization problem leads to a family of linear equations. For example, the maximum flow and 
minimum cost flow problems require the solution of linear systems whose matrices are symmetric 
and diagonally-dominant. The generalized versions of these flow problems result in symmetric 
M- matrices. 

The generalized maximum flow problem is specified by a directed graph (V, E), an inward 
capacity c(e) > and a multiplier 7(e) > for each edge e, and source and sink vertices s and t. 
For every unit flowing into edge e, 7(e) flows out. In lossy generalized flow problems, each multiplier 
7(e) is restricted to be at most 1. In the generalized maximum flow problem, one is asked to find 
the flow / : E — > H + that maximizes the flow into t given an unlimited supply at s, subject to 
the capacity constraints on the amount of flow entering each edge. In the generalized minimum 
cost flow problem, one also has a cost function q(e) > 0, and is asked to find the maximum flow of 
minimum cost (see [AM093] ). 

In the following chart, we compare the complexity of our algorithms with the fastest algorithms 
of which we are aware. The running times are given for networks in which all capacities and costs 
are positive integers less than U and every loss factor is a ratio of two integers less than U. For 
the standard flow problems, our algorithms are exact, but for the generalized flow problems our 
algorithms find additive e approximations, while the other approximation algorithms have multi- 
plicative error (1 + e). However, we note that our algorithms only require arithmetic with numbers 
of bit-length O (\og(nU/e)), whereas we suspect that the algorithms obtaining multiplicative ap- 
proximations might require much longer numbers. 



In the chart, C refers to the value of the flow. 



Exact algorithms 


Approximation algorithms 


Our algorithm 


Generalized Maximum Flow 

O (m 2 (m + n log n) log U) [G.TQ97] 
O (m L5 n 2 log(ncO) [Vai89| 


O (m 2 /e 2 ) [FW02] 

O (m{m + n log log B) log e _1 ) 

[GFNR98] [TW98J |FW02| 


O (m l - 5 log 2 (U/e)) 


Generalized Minimum Cost Flow 

O (m l - b n 2 \og{nU)) |Vai89j 


0(m 2 loglogB/e 2 ) [FW02J 


O {m 1 - 5 'log 2 '(U/e)) 


Maximum Flow 

O (min(n 3//2 , m 1 / 2 )mlog(n 2 /m) log U) 
|GR98j 




6 (m L5 log 2 C/) 


Minimum Cost Flow 

O {nm\og{n 2 /m)\og(nC)) fGT87j 
O (nmtioe Joe U) \oe(nC)) [AGOT92] 
O ((m log n)(m + n log n)) [Orl88] 




6 (m L5 log 2 C/) 



1.1 The solution of systems in M-matrices 

A symmetric matrix M is diagonally dominant if each diagonal is at least the sum of the absolute 
values of the other entries in its row. A symmetric matrix M is an M -matrix if there is a positive di- 
agonal matrix D for which DMD is diagonally dominant. Spielman and Teng [ST04, ST06] showed 
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how to solve linear systems in diagonally dominant matrices to e accuracy in time O (mloge -1 ). 
We show how to solve linear systems in M-matrices by first computing a diagonal matrix D for 
which DMD is diagonally dominant, and then applying the solver of Spielman and Teng. Our 
algorithm for finding the matrix D applies the solver of Spielman and Teng an expected O (log n) 
times. While iterative algorithms are known that eventually produce such a diagonal matrix D, 
they have no satisfactory complexity analysis jLi02|. |LLH + 98~1 IBCPT05] . 

1.2 Analysis of interior point methods 

In our analysis of interior-point methods, we examine the complexity of the short-step dual path 
following algorithm of Renegar [Ren88j as analyzed by Ye [Ye97| . The key observations required by 
our complexity analysis are that none of the slack variables become too small during the course of 
the algorithm and that the algorithm still works if one O (l/-y/m)-approximately solves each linear 
system in the matrix norm (defined below). Conveniently, this is the same type of approximation 
produced by our algorithm and that of Spielman and Teng. This is a very crude level of approxi- 
mation, and it means that these algorithms can be applied very quickly. While other analyses of 
the behavior of interior point methods with inexact solvers have appeared |Ren96j , we are unaware 
of any analyses that are sufficiently fine for our purposes. 
This analysis is given in detail in Appendix [Cl 

1.3 Outline of the paper 

In Section [21 we describe the results of our analysis of interior point methods using apprximate 
solvers. In Section [3l we describe the formulation of the generalized flow problems as linear pro- 
grams, and discuss how to obtain the solutions from the output of an interior-point algorithm. In 
Section^ we give our algorithm for solving linear systems in M-matrices. 

2 Interior-Point Algorithm using an Approximate Solver 

Our algorithm uses numerical methods to solve a linear program formulation of the generalized flow 
problems. The fastest interior-point methods for linear programs, such as that of Renegar |Ren88| 
require only O (\/n) iterations to approach the solution, where each iteration takes a step through 
the convex polytope by solving a system of linear equations. 

In this paper, we consider stepping through the linear program using an only an approximate 
solver, i.e. an algorithm x = Solve(M, 6,e) that returns a solution satisfying 

\\x-M-H\\ M <e\\M^b\\ M 

where the matrix norm is given by = V v T Mv. 

As mentioned above, we have analyzed the Renegar {Ren88j version of the dual path-folllowing 
algorithm, along the lines of the analysis that found in |Ye97j , but modified to account for the use 
of an approximate solver. 

In particular, using the approximate solver we implement an interior-point algorithm with the 
following properties: 

Theorem 2.1. x = InteriorPoint(^4, 6, c, A m j n , T, e) takes input that satisfies 

• A is an n x m matrix; 

b is a length n vector; c is a length m vector 
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• AA T is positive definite, and X m in > is a lower bound on the eigenvalues of AA T 

• T > is an upper bound on the absolute values of the coordinates in the dual linear program, 
i.e. 

WyW^ <T and IML^T 

for all (y, s) that satisfy s = c — A T y > 

• initial point y° is a length n vector where A T y° < c 

• error parameter e satisfies < e < 1 

and returns x > such that \\Ax — b\\ < e and c T x < z* + e. 
Let us define 

• U is the largest absolute value of any entry in A,b, c 

• s^ in is the smallest entry of s° = c — A T y° 

Then the algorithm makes O ( ^pm log Tu jp- — ) calls to the approximate solver, of the form 

Solve (AS~ 2 A T + vv T , ■, e') 

where S is a positive diagonal matrix with condition number O u e m ' ^ , and v,e' satisfy 

H| / TUm\ 

log = 0M g 

c V *mm fc / 

In Appendix £3 we present a complete description of this algorithm, with analysis and proof of 
correctness. 



3 Solving Generalized Flow 

We consider network flows on a directed graph (V,E) with V = [n], E = {ei,--- ,e m }, source 
s £ V and sink t £ V. Edge goes from vertex Vj to vertex Wj. and has inward capacity c(ej), 
flow multiplier j(ej) < 1, and cost q(ej). 

We assume without loss of generality that t has a single in-edge, which we denote as et, and no 
out-edges. 

The generalized max-flow approximation algorithm will produce a flow that sends no worse 
than e less than the maximum possible flow to the sink. 

The generalized min-cost approximation algorithm will produce a flow that, in addition to being 
within e of a maximum flow, also has cost no greater than the minimum cost of a maximum flow 
(see [FW02] ). 



3.1 Fixing Approximate Flows 

The interior-point algorithm described in the previous section produces an output that may not 
exactly satisfy the linear constraints Ax = b. In particular, when we apply the algorithm to a 
network flow linear program, the output may only be an approximate flow: 

Definition 3.1. An e-approximate flow approximately satisfies all capacity constraints and flow 
conservation constraints. In particular, every edge may have flow up to e over capacity, and every 
vertex besides s and t may have up to e excess or deficit flow. 

An exact flow satisfies all capacity constraints and has exact flow conservation at all vertices 
except s and t. 
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We are going to modify the graph slightly before running the interior-point algorithm, so that 
it will be easier to obtain an exact flow from the approximate flow given by the interior-point 
algorithm. 

Let us compute the least-lossy-paths tree T rooted at s. This is the tree that contains, for each 
v £ V — {s,t}, the path tt s ^ v from s to d that minimizes L(v) = Y\ e€lTs „ 7( e ) _1 ) the factor by 

which the flow along the path is diminished. We can find this tree in time O (m), using Dijkstra's 
algorithm to solve the single-source shortest-paths problem with edge weights — log7(e). 

Next, we delete from the graph all vertices v such that L{v) > 2 m n u ■ Note that in a maximum- 
flow, it is not possible to have more than ^ flowing into such a v, since at most mil can flow out 
of s. Thus, deleting each such v cannot decrease the value of the maximum flow by more than 
In total, we may decrease the value of the maximum flow by at most | . 

2 

We define cflow = ^^i^u'^ ' ^ n ^ e subsequent sections, we show how to use the interior-point 
method to obtain an e^cw-approximate how that has a value within €4 of the maximum flow. 
Assuming that the graph had been preprocessed as above, we may convert the approximate flow 
into an exact flow: 

Lemma 3.2. Suppose all vertices v £ V — {s,t} satisfy L(v) < 2 mnU • ^ n ® ( m ) ^ me > we are a ^ e 
to convert an €FLOW- a PP rox 'i' ma t e flow that has a value within | of the maximum flow into an exact 
flow that has a value within | of the maximum flow. The cost of this exact flow is no greater than 
the cost of the approximate flow. 

Proof. Let us first fix the flows so that no vertex has more flow out than in. We use the least-lossy- 
paths tree T, starting at the leaves of the tree and working towards s. To balance the flow at a 
vertex v we increase the flow on the tree edge into v. After completing this process, for each v we 

2 

will have added a path of flow that delivers at most Q 4m ^ n -2 U -i additional units of flow to v. Since 
L(v)< 

2mnU > no suc h path requires more than g4 m 2 n 2(/3 ' € ^ — 32rnnU 2 ^ ow on an e< ^S e ' an d so 
in total we have added no more than 32 mU' 2 ^° eacn edge. 

Next, let us fix the flows so that no vertex has more flow in than out. We follow a similar 
procedure as above, except now we may use any spanning tree rooted at and directed towards t. 
Starting from the leaves, we balance the vertices by increasing flow out the tree edge. Since the 

2 2 

network is lossy, the total amount added to each edge is at most 64m f w2c/ 3 ■ n < Q 4n ^i nU i ■ 

2 

Recall that we started with each edge having flow up to 64Tn 2 w 2;y3 over capacity. After balancing 
the flows at the vertices, each edge may now be over capacity by as much as 

e 2 e 2 



32mU 2 6Am 2 nU 3 64m 2 n 2 U 3 ~ IQmU 2 

Since the edge capacities are at least 1, the flow on an edge may be as much as (1 + j^fi) times 
the capacity. 

Furthermore, while balancing the flows we may have added as much as 16r ^ 2 ■ thU = -^j to 
the total cost of the flow. Assuming that the value of approximate flow was at least |, its cost 
must also have been at least f, and so we have increased the cost by a multiplicative factor of at 
most (1 + ^j). 

(If the approximate flow had value less than |, then the empty flow trivially solves this flow 
rounding problem.) 

By scaling the entire flow down by a multiplicative factor of (1 + gfj) -1 , we solve the capacity 
violations, and also reduce the cost of the exact flow to be no greater than that of the approximate 
flow. Since the value of a flow can be at most U, the flow scaling decreases the value of the flow by 
no more than e/4, as required. □ 
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The above procedure produces an exact flow that is within e/2 of the maximum flow in the 
preprocessed graph, and therefore is within e of the maximum flow in the original graph. Further- 
more, the cost of the flow is no greater than the minimum cost of a maximum flow in the original 
graph. 

Thus to solve a generalized flow problem, it remains for us to describe how to use the interior- 
point algorithm to generate a ei^o^- a PP rox i ma te flow that has a value within e/4 of the maximum 
flow, and, for the min-cost problem, also has cost no greater than the the minimum cost of a 
maximum flow. 



3.2 Generalized Max-Flow 

We formulate the maximum flow problem as a linear program as follows: Let A be the (n — 2) x m 
matrix whose nonzero entries are A v .j = —1 and A w . j = j(ej), but without rows corrsponding 
to s and t. Let c be the length m vector containing the edge capcities. Let u t be the length m 
unit vector with a 1 entry for edge ej. Let the vectors x\ and X2 respectively denote the flow into 
each edge and the unused inward capacity of each edge. The max-flow linear program, in canonical 
form, is: 



mm — u t x\ 

Xi 



S.t. 



'A 




X\ 




"0" 


I I 




. X 2. 




c 



and Xi > 



The constraint Ax\ = ensures that flow is conserved at every vertex except s and t, while the 
constraint x\ + X2 = c ensures that the capacities are obeyed. 

Now, the dual of the above linear program is not bounded, which is a problem for our interior- 
point algorithm. To fix this, we modify the linear program slightly: 



• , t 4 ^ 
mm — u t Xi H 

x i V tFLOW 



(1^X 3 + ln_2 x i + 1 n-2 :E 5) 



and 



A I 
I I -I 



s.t 


X 


i > 




~x{ 






I 




X2 




"0" 






X3 














c 






X4 








X5. 







(We use lfc to denote the all-ones vector of length k.) 

Lemma 3.3. This modified linear program has the same optimum value as the original linear 
program. 

Proof. Let us examine the new variables in the modified program and note that X3 has the effect 
of modifying the capacities, while £4 and x§ create excess or deficit of flow at the vertices. Since 
we have a lossy network, a unit modification of any of these values cannot change the value of the 
flow by more than 1, and therefore must increase the value of the modified linear program. Thus, 
at the optimum we have X3 = X4 = x^ = and so the solution is the same as that of the original 
linear program. □ 
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The modified linear program has the following equivalent dual linear program: 

s.t. Sj > 



T 

max c y 2 



and 





I ' 




~Sl 






I 








S2 






-I 




+ 


S3 




I 










S 4 




—I 






, S 5_ 







(4:11 /e F low) ■ J -m 
{AU/e FL ow) ■ ln-2 



Lemma 3.4. The above dual linear program is bounded. In particular, the coordinates of all feasible 



dual points have absolute value at most (nil + 1) • 



tFLOW 



+ 1. 



Proof. Of the five constraints in the dual linear program, the last four give 



W 



as an explicit 



CFLOW 

bound on the absolute value of y coordinates. It then follows that — — — is a upper bound on 



the coordinates of 82, £3,84,85, and the coordinates of Si 



(nU + 1) 



W 



CFLOW 



+ 1. 



y 2 can be at most 
□ 



We refer to the Sj variables as the slacks. Recall that we must provide the interior-point 
algorithm with an initial dual feasible point y° such that the corresponding slacks s° are bounded 
away from zero. We choose the following initial point, and note that the slacks are bounded from 
below by 



tFLOW 
















_—(2U/eFLOw) ■ l?n_ 








(2U /eFhow) ■ 1m _ 
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[2U /tFLOw) ■ 1m 


,0 




(2U/eFLOw) ■ 1m 


,0 




{4U/e F LOw) ■ ln- 


-2 


S° 
- s 5. 




(4U/eflOw) ■ ln- 


-2 



We must also provide the interior-point algorithm with a lower bound on the eigenvalues of the 
matrix 

' A T I 
I 

-I 

I 

-I 



A I -I 

I I -I 



AA T + 2I A' 
A T 31 



Note that we may subtract 21 from the above matrix and still have a positive definite matrix, so 
Amm = 2 is certainly a lower bound on the eigenvalues. 

Using the above values for y° and \ m in, an d the bound on the dual coordinates given in Lemma 
13.41 we now call InteriorPoint on the modified max-flow linear program, using error parameter 
zflow _ j n solution returned by the interior-point algorithm, the vector x\ assigns a flow value 
to each edge such that the flow constraints are nearly satisfied: 



Lemma 3.5. x\ is an eFLOW- a PP rox ^ ma t e fl ow with value within €flow/2 of the maximum flow. 
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Proof. Observe that the amount flowing into t is at least —1 times the value of the modified linear 
program. Since the interior-point algorithm generates a solution to the modified linear program 
within eFLOw/2 of the optimum value, which is —1 times the maximum flow, the amount flowing 
into t surely must be within e flow I '2 of the maximum flow. 

Now, let us note more precisely that the modified linear program aims to minimize the objective 
function computed by subtracting the amount flowing into t from 4U INFLOW times the sum of the 
entries of 333, 334, and 335. Since the minimum value of this objective function must be negative, and 
the solution returned by the interior-point algorithm has a value within €flow/2 of the minimum, 
the value of this solution must be less than €flow/2 < U. The amount flowing into t is also at 
most U, so no entry of 333,334,335 can be greater than 2U / \AU / 'cflow) = cflow/2- 

The interior-point algorithm guarantees that 



H-Axi + 334 — 335 1| < 
and so we may conclude that 



(■FLOW 



and 



\X\ + 332 — 333 — C < 



(FLOW 



\\Ax\\\ < cflow and x\ < c + cflow 

Indeed, this is precisely what is means for x± to describe an e^oty-approximate flow. 



□ 



3.3 Generalized Min-Cost Flow 

As a first step in solving the generlized min-cost flow problem, we solve the generalized max-flow 
linear program as described above, to find a value F that is within | of the maximum flow. 

We now formulate a linear program for finding the minimum cost flow that delivers F units of 
flow to t: 

T 

mm q X\ 

Xi 

and 





s.t. 


33 


i > 


A 






33i 




F ■ e t ~ 


I 


I 




X 2_ 




c 



where q is the length n vector containing the edge costs, and e< is the length n — 1 vector that 
assigns 1 to vertex t and to all the other vertices except s. A is the same matrix as in the 
max-flow linear program, except that we include the row corresponding to t, which translates to a 
new constraint that F units must flow into t. 

We must again modify the linear program so that the dual will be bounded: 



T ( 4mU 2 
min I q x\ + 

xi \ \£FLOW 



) (l^3 + lLl*4 + lLl*5)) 



and 



A I 
I I -I 



s.t 


33 


i > 




~x{ 






I 




x 2 




F ■ ei 






X 3 












c 






334 








_33 5 _ 







Lemma 3.6. This modified linear program has the same optimum value as the original linear 
program. 
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Proof. We examine the new variables and note that 213 modifies the capacities, while 2:4 and x§ 
create excess supply (or demand) at the vertices. A unit modification to any of these values can at 
best create a new path for one unit of flow to arrive at the sink. This new path has cost at least 1, 
and it can replace an path in the optimum flow of cost at most nU, for a net improvement in the 
cost of the flow of at most nU — 1, which is less than AmU . Thus the value of the modified linear 

(FLOW 

program can only increase when these new variables are set to non-zero values. □ 



Now, the dual linear program is: 



max (F ■ efy 1 + c T y 2 ) 



and 



'A 1 



I ' 






I 




Vi 






-I 




+ 


S3 






.2/2. 




s 4 












S5 



s.t. Si > 

q 


(4m£7 2 Inflow) ■ l m 
{AmU 2 /e FLO w) ■ ln-l 
(4mU 2 /e FL ow) ■ ln-i. 



Lemma 3.7. The above dual linear program is bounded. In particular, the coordinates of all feasible 
dual points have absolute value at most (nU + 1) • imU . 

' (FLOW 

program, the last foi 
bound on the absolute value of y coordinates. It then follows that N ' " 



(FLOW 

Proof. Of the five constraints in the dual linear program, the last four give 

r2 



(FLOW 



4mU as an explicit 

(FLOW 

is a upper bound on 



the coordinates of 82, S3, 84, S5, and the coordinates of s\ = q — A T y 1 — y 2 can be at most 
(nU + 1) 



4mU 2 
(-FLOW ' 

Let us also note that y° 

mil 2 



-(mU 2 /e FL ow)h 



□ 



is an initial interior dual point with all slacks 



^FLOW 



at least 

Using the above initial point, the bound on the dual coordinates from Lemma [3.7l and \ r 
as in the previous section, we run InteriorPoint on the modified min-cost linear program, with 
error parameter eph ° w . In the solution returned by the interior-point algorithm, the vector x\ 
assigns a flow value to each edge such that the flow constraints are nearly satisfied: 

Lemma 3.8. x\ is an eF F ow- a VV T0X ^ ma ^ e fl ow with value within || of the maximum flow. 

Proof. Note that any flow in total cannot cost more that mil 2 , even if all edges are filled to 
maximum capacity. Therefore the value of the solution output by the interior-point algorithm can 
be at most mil 2 + €FL 2 olv < 2mU 2 , and so in particular no entry of £3, 214, x§ can be greater than 

(FLOW 
2 



Now, the interior-point algorithm guarantees that 

cflow 



\\Ax\ + ;e 4 — £5 — F ■ et\\ < 
and so we may conclude that 

\\Axi - F ■ e t \\ < (-Flow 



and 



and 



35 1 + x 2 - S3 



c < 



tFLOW 



Xi < C + 6FLOW 



These inequalities imply that this is a €flow- approximate flow, and additionally that at least 
F — cflow is flowing into t. Since F is within I of the maximum flow, the amount flowing into t 



must be within f + €flow < f§ of the maximum flow. 



□ 
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By scaling down the x\ flow slightly, we obtain a flow that does not exceed the minimum cost 
of a maximum flow: 

Lemma 3.9. x' x = (1 — -^j)x\ is an epi J ow- a VV TOX ' l ' ma ^ e fl ow with value within | of the maximum 
flow, and with cost at most the minimum cost of a maximum flow. 

Proof. We may assume that the value of flow x\ is at least |f , because otherwise the maximum 
flow would have to be at most || + || = | , and so the empty flow would trivially be within | of 
the maximum. Therefore, the minimum cost of a maximum flow must also at be least |f . 

The interior-point algorithm guarantees that the cost of X\ does not exceed this optimum cost 
by more than eFL ° w ; and so must also not exceed the optimum cost by a multiplicative factor of 
more than (1 + 16ej ^ cw ) < (1 + Thus. x' x = (1 — j§jj)xi must have cost below the optimum. 

Furthermore, since the value of the flow x± can be at most U, scaling down by (1 — j^jj) cannot 
decrease the value of the flow by more than j^. Therefore, the value of the value x± is within 
■j2 + §§ < f of the maximum. □ 



3.4 Running Time 

The linear systems in the above linear programs take the form 



A 



A I -I 

I I -I 



so the running time of the interior-point method depends on our ability to approximately solve 
systems of the form AS~ 2 A T + vv T , where diagonal matrix S and vector v are as described in 
Theorem 12.11 As it turns out, this is not much more difficult than solving a linear system in 
AS 'i A T , where S± is the upper left submatrix of S. 

The matrix AS^ 2 A T is a symmetric M-matrix. In the next section, we describe how to ap- 
proximately solve systems in such matrices in expected time O (mlog where k is the condition 
number of the matrix. We then extend this result to solve the systems AS~ 2 A T + vv T in time 
O (mlog ] , where k is the condition number of AS~ 2 A T . 



Theorem 3.10. Using out interior-point algorithm, we can solve the generalized max flow and 
generalized min-cost flow problems in time O (m 3 / 2 log 2 ([//e)) 

Proof. According to Theorem l2.14 the interior-point algorithm requires O ( \/mlog , ^Jy 1 — ) calls 
to the solver. 

Recall that T is an bound on the coordinates of the dual linear program, and is the smallest 
slack at the initial point. Above, we gave both of these values to be polynomial in for both the 
max-flow and min-cost linear programs. We also gave \ m in = 2 as a lower bound on the eigenvalues 
of AA T . Thus, the total number of solves is O (\/m log — ) . 

Again referring to Theorem 12. 1\ we find that the condition number of AS~ 2 A T is be polynomial 
in as is the expression We conclude that each solve takes time O (mlog — ). 

The preprocessing only took time O (m) so we obtain a total running time of O (m 3 / 2 log 2 (?7/e)) . 

□ 



3.5 Standard Min-Cost Flow 

In this section we describe how to use interior-point algorithms to give an exact solution to the 
standard (i.e. no multipliers on edges) min-cost flow problem. 
We use the following property of the standard flow problem: 
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Theorem 3.11 (see [Sch03( Theorem 13.20]). Given a flow network with integer capacities, and a 
positive integer F , let SIflow be the set of flow vectors x that flow F units into t and satisfy all 
capacity and flow conservation constraints. Then Sip LOW i> s a convex poly tope in which all vertices 
have integer coordinates. 

Our goal is to find the flow in Sip low of minimum cost. Since the cost function is linear, if 
there is a unique minimum-cost flow of value F, it must occur at a vertex of SI flow- By Theorem 
13. Ill this must be an integer flow, and we could find this flow exactly by running the interior-point 
algorithm until it is clear to which integer flow we are converging. 

Unfortunately, the minimum-cost flow may not be unique. However, by applying the Isolation 
Lemma of Mulmuley, Vazirani, and Vazarani [MVV87] , we can modify the cost function slightly 
so that the minimum-cost flow is unique, and is also a minimum-cost flow under the original cost 
function. 

Let us first state a modified version of the Isolation Lemma: 

Lemma 3.12 (sec [KS01, Lemma 4]). Given any collection of linear functions on m variables 
with integer cooefficients in the range {0, . . . , U}. If each variable is independently set uniformly at 
random to a value from the set {0, . . . ,2mU}, then with probability at least 1/2 there is a unique 
function in the collection that takes minumum value. 

We now describe how to force the minimum-cost flow to be unique: 

Lemma 3.13. Given a flow network with capacities and costs in the set {1, 2, . . . , U}, and a positive 
integer F, modify the cost of each edge independently by adding a number uniformly at random from 
the set { irj^u' 2 > 4n$U' 2 ' ' ' ' ' Im^ZT 7 }' ^ en with probability at least 1/2, the modified network has a 
unique minimum-cost flow of value F, and this flow is also a minimum-cost flow of value F in the 
original network. 

Proof. The modified cost of a flow at a vertex of Sip low is a linear function of m independent vari- 
ables chosen uniformly at random from the set { , zJtjj'i > • • • > ^S\j 2 } • where the coefficients 
are the coordinates of the flow vector, which by Lemma [3.11l are integers in the range {0, . . . , U}. So 
the Isolation Lemma tells us that with probablity at least 1/2, there is a unique vertex of SI flow 
with minumum modified cost. 

Now, any vertex that was not originally of minimum cost must have been more expensive than 
the minimum cost by an integer. Since the sum of the flows on all edges can be at most mU, and 
no edge had its cost increased by more than , the total cost of any flow cannot have increased 
by more than 1/2. Thus, a vertex that was not originally of minimum cost cannot have minimum 
modified cost. □ 

We may now give an exact algorithm for standard minimum-cost flow. Note that this algorithm 
works for any integer flow value, but in particular we may easily find the exact max-flow value by 
running the interior-point max-flow algorithm with an error of 1/2, since we know the max-flow 
value is an integer. 

Lemma 3.14. To solve the standard minimum-cost flow problem in expected time O (m 3//2 log 2 U) , 
perturb the edge costs as in Lemma \3.13l then run the min- cost flow interior point algorithm with 
an error of 12 m a £/ a , and round the flow on each edge to the nearest integer. 

Proof. Let us prove correctness assuming that the modified costs do isolate a unique minimum-cost 
flow. The running time then follows directly from Theorem l3.im and the fact from Lemma r3.13l that 
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after a constant number of tries we can expect the modified costs to yield a unique minimum-cost 
flow. 

We first note that the modified edge costs are integer multiples of S = ^rfp ■ Therefore, by 
Theorem 13,111 the cost of the minumum-cost flow is at least 5 less than the cost at any other vertex 
of £Iflow- 

Now, the flow returned by the interior-point algorithm can be expressed as a weighted average 
of the vertices of Qflow- Since the cost of this flow is within i 2rn ^t/3 = gfj of the minimum 
cost, this weighted average must assign a combined weight of at most to the non- minimum-cost 
vertices. Therefore, the flow along any edge differs by at most 1/3 from the minimum-cost flow. 
So by rounding to the nearest integer flow, we obtain the minimum-cost flow. □ 



4 Solving linear systems in symmetric M-Matrices 

A symmetric M -matrix is a positive definite symmetric matrix with non-positive off-diagonals (see, 
e.g. [HJ91, Axe96, BP94J). Every M-matrix has a factorization of the form M = AA T where each 
column of A has at most 2 nonzero entries [BCPT05 . Given such a factorization of an M-matrix, 
we we will show how to solve linear systems in the M-matrix in nearly-linear time. Throughout 
this section, M will be an n x n symmetric M-matrix and A will be a n x m matrix with 2 nonzero 
entries per column such that M = AA T . Note that M has O (m) non-zero entries. 

Our algorithm will make use of the Spielman-Teng O (m) expected time approximate solver 
for linear systems in symmetric diagonally-dominant matrices, where we recall that a symmetric 
matrix is diagonally-dominant if each diagonal is at least the sum of the absolute values of the other 
entries in its row. It is strictly diagonally- dominant if each diagonal execceds each corresponding 
sum. 

We will use the following standard facts about symmetric M-matrices, which can be found, for 
example, in [HJ91]: 



Fact 4.1. If M 



Mn M12 
Mf 2 M 22 



is a symmetric M-matrix with Mn a principal minor, then: 



1. M is invertible and M~ 1 is a nonnegative matrix. 

2. Mi 2 is a nonpositive matrix. 

3. Mn is an M-matrix. 

4. The Schur complement S = M 22 - M^M^Mu is an M -matrix. 

5. If all eigenvalues of M fall in the range [A m ^, A m aa:]? then so do all diagonal entries of S. 

6. For any positive diagonal matrix D, DMD is an M-matrix. 

7. There exists a positive diagonal matrix D such that DMD is strictly diagonally-dominant. 

Our algorithm will work by finding a diagonal matrix D for which DMD is diagonally-dominant, 
providing us with a system to which we may apply the solver of Spielman and Teng. Our algorithm 
builds D by an iterative process. In each iteration, it decreases the number of rows that are 
not dominated by their diagonals by an expected constant factor. The main step of each iteration 
involves the solution of O (log n) diagonally-dominant linear systems. For simplicity, we first explain 
how our algorithm would work if we made use of an algorithm x = ExactSolve(M, b) that exactly 
solves the system Mx = b, for diagonally-dominant M. We then explain how we may substitute 
an approximate solver. 

The key to our analysis is the following lemma, which says that if we multiply an M-matrix by 
a random diagonal matrix, then a constant fraction of the diagonals probably dominate their rows. 
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Lemma 4.2 (Random Scaling Lemma). Given an n x n M -matrix M, and positive real values 
( < 1 and r < j, let D be a random diagonal n x n matrix where each diagonal entry d{ is chosen 
independently and uniformly from the interval (0, 1). 

Let T C [n] be the set of rows of MD with sums at least r times the pre-scaled diagonal, i.e. 

T = {% g [n] : (MDl)j > rmu} 

With probability at least \~ty , we have 



1 r 

1 



1 



n 



where (3 is the fraction of the diagonal entries of M that are less than ( times the average diagonal 
entry. 

Note in particular that for r = 0, T is the set of rows dominated by their diagonals. 
We will use the Random Scaling Lemma to decrease the number of rows that are not dominated 
by their diagonals. We will do this by preserving the rows that are dominated by their diagonals, 

"Mn M 12 



and applying this lemma to the rest. Without loss of generality we write M 



A 2 A\ A 2 Al 



M 12 



M 22 



where the rows in the top section of M are the ones that are already diagonally- 



dominant, so in particular Mn is diagonally-dominant. Let S = M 22 — Mj^Mj^ M\ 2 be the Schur 
complement and let Sd be the matrix containing only the diagonal entries of S. 

We construct a random diagonal matrix Dr of the same size as M 22 by choosing each diagonal 

ID! 



element independently and uniformly from (0, 1). We then create diagonal matrix D 



-M u L M l2 D 2 l. 



D 2 

We know that 



— 1/2 

where D 2 = S D Dr and the diagonal entries of D\ are given by 

the diagonal entries of D\ are positive because Fact 14.11 tells us that MQ 1 is nonnegative and M12 
is nonpositive. 

We now show that the first set of rows of DMD are diagonally-dominant, and a constant fraction 
of the rest probably become so as well. Since M is an M-matrix and D is positive diagonal, DMD 
has no positive off-diagonals. Therefore, the diagonally-dominant rows of DMD are the rows with 
nonnegative row sums. The row sums of DMD are: 



DMD1 



D 1 M U D 1 1 + D X M 12 D 2 1 
D 2 Mf 2 D x l + D 2 M 22 D 2 \ 





D 2 SD 2 1 





D R S~ 1/2 SS n 1/2 D R l 



Note that the diagonal entries of Sp 1 ^ 2 SS D 1/Z are all 1. Thus by invoking Lemma 14.21 with r = 
and (" = 1, we find that there is a 1/7 probability that at least 1/24 of the row sums in the bottom 
section of DMD become nonnegative. Furthermore, we see that row sums in the top section remain 
nonnegative. 

The only problem with this idea is that in each iteration it could take O (mn) time to compute 
the entire matrix S. Fortunately, we actually only need to compute the diagonals of S, (i.e. the 
matrix Sd). In fact, we only actually need a diagonal matrix S that approximates Sd- As long 
the diagonals of T,~ l / 2 SH~ 1 / 2 fall in a relatively narrow range, we can still use the Random Scaling 
Lemma to get a constant fraction of improvement at each iteration. 



-1/2 
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To compute these approximate diagonal values quickly, we use the random projection technique 
of Johnson and Lindenstrauss |JL84j . In AppendixEJ we prove the following variant of their result, 
that deals with random projections into a space of constant dimension: 

Theorem 4.3. For all constants a,(3,^,p E (0,1), there is a positive integer k = kjL (en, /3, 7, p) 
such that the following holds: 

For any vectors V\, ... ,v n €. M m let R be a k x m matrix with entries chosen independently at 

random from the standard normal distribution, and let W{ = ^J^Rvi. 
With probability at least p both of the following hold: 



I- <(l + 7)n 



2. 



< fin 



II- 

Let us note that 

S = A 2 {I - AjM^A^A^ = A 2 (I - A\M^A x fA\ 

because (/ — Af M^ 1 A\) is a projection matrix. So if we let ctj denote the ith row of A 2 , we can 
write the ith diagonal of S as su = ||(J — AjM^A^afW 2 . Then if we use Theorem 14.31 to create a 
random projection matrix R, \\R(I — A\M^A\)aJ\\ 2 gives a good approximation to sa. Moreover, 
we can use one call to ExactSolve to compute each of the constant number of rows of the matrix 
P = R(I — AjM^ A\). Since A 2 has O (m) entries, we can compute PA\ in O (m) time, and 
obtain all the approximations ||Pa^|| 2 in O (m) time, yielding the desired approximations of all sa 
values. 

Our suggested algorithm, still using an exact solver, is given in Figure [H To make this algorithm 
fast, we replace the calls to the exact solver with calls to the approximate solver STSolve of Spielman 
and Teng: 

Theorem 4.4 (Spielman- Teng [ST04} IST06] ). The algorithm x = STSolve(M, b,e) takes as input 
a symmetric diagonally- dominant n x n matrix M with m non-zeros, a column vector b, and an 
error parameter e > 0, and returns in expected time O (mlog(l/e)) a column vector x satisfying 

\\x - M _1 6|| A/ < e||M _1 6|| M 

We define the algorithm MMatrixSolve(A, 6, e, A m j n , \ m ax) as a modification of the algorithm 
ExactMatrixSolve in Figure [H For this algorithm we need to provide upper and lower bounds 
A mflI , Xmin on the eigenvalues of the matrix A, and the running time will depend on k = Xmax/^min- 

The modifications are that we need to set parameters: 



5 = (l/2A)\]£ n K- 1/2 n- 1 ei = .005(1.0lKmn)- 1 / 2 e 2 = (l/72)^ 5 / 2 n 
and substitute the calls to ExactSolve in lines 2c, 2h and 3 respectively with 

• S1Solve(D 1 M 11 D 1 ,D 1 A 1 rf,e 1 ) 

• STSolvepiMnDi, -Di(-Mi2£>2 + e 2 ) 

• S7Solve(DMD,Db,e). 

We may note that the final call to STSolve guarantees that 

\\D- l x - D- l M- l b\\ DM D < e\\D- 1 M- 1 b\\ D MD 
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x = ExactMMatrixSolve(A, b) 
Given: n x m matrix A, where M 
column. 

Returns: x satisfying Mx = b 



AA T is an M-matrix and A has at most 2 non-zeros per 



1. Set D := I. 

2. Until DMD is diagonally dominant do: 

a. Permute so that DMD - ^ ^' 1 1 ^ 1 ^' Jii ' 



D 1 A 1 A\D 1 D x A x A\Di 
D 2 A 2 A\D X D 2 A 2 A\D 2 



has 



D 2 Mj 2 D l D 2 M 22 D 2 
the diagonally dominant rows in the top section. Let &\i . . . o>ii be the rows of A 2 . 

b. Set k = kjL(ijjQ, I' TSo' i)' an( ^ ^ be a random k x m matrix with independent 
standard normal entries. Let be the ith row of R. 

c. For i = 1, . . . , k, compute qf = ExactSolve(L>iMnZ?i, D\A\rJ). 

T\\2 



d. Set Q= [qj ■■■ q+ kj 

e. Let S be the v x v diagonal matrix with entries Oi 



\\{R-QDxAi)ai 

f. Let Dr be a uniform random v x v diagonal matrix with diagonal entries in (0, 1). 

g. Set D' 2 = TrV*D R 

h. Set D[ to be the matrix with diagonal D\ ■ ExactSolve(Z?iMiiDi, — DiMi 2 D' 2 l) 



Set D 



D' 



3. Return x = D ■ ExactSolve(DMD, Db) 



Figure 1: Algorithm for solving a linear system in a symmetric M-matrix. To speed up the algorithm 
we will replace the exact solver with the Spielman Teng approximate solver. 



or equivalently 



x 



M~ 1 b\\ M < ellM^bl 



M 



so the output fulfills the specification of an approximate solver, provided that the algorithm termi- 
nates. 

We can in fact bound the running time of this algorithm as follows: 

Theorem 4.5. The expected running time of the algorithm MMatrixSolve is O (mlog^). 

Proof. The running time is dominated by the calls to the Spielman- Teng solver. There are O (1) 
such solves per iterations, each of which take time O (mlog«;), and at the conclusion of the algo- 
rithm, there is one final call of time O (mloge -1 ). 

So, to prove the running time, it suffices for us to give a O (log m) bound on the expected 
number of iterations. In particular, it suffices to show that in each iteration, the number of non- 
diagonally-domainant rows in DMD decreases by a constant fraction with constant probability. 

\D l 

D 2 



In analyzing a single iteration, we let D 



D' 



of the iteration, and we let D' 
prove: 

Lemma 4.6. D' is a positive diagonal matrix. 



denote the diagonal scaling at the start 
denote the new diagonal scaling. In Appendix |Aj we 
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This implies that D'MD' has no positive off-diagonals, thereby enabling us to check which rows 
of D'MD' are diagonally-dominant by looking for rows with nonnegative row sums. 

We again let S = M22 — Mj^M-Q 1 M12 denote the Schur complement, and let So denote the 
matrix containing the diagonal entries of S. Let us also define S = E -1 / 2 ^!]" 1 / 2 . We know from 
Facts ETT14 and[3j]6 that § is an M-matrix. 

Let Sd be the matrix containing the diagonal entries of S. In Appendix [Aj we show that the 
row sums of MD' are related to S as follows: 

Lemma 4.7. 

md,1 -[tM\sd r -\~s d )i 

The upper part of the above inequality tells us that all the row sums that were nonnegative in 
DMD remain nonnegative in D'MD' . From the lower part of the inequality and by invoking the 
Random Scaling Lemma on the matrix S with r = g, we find that with probabilty at least ^3, the 

fraction of remaining rows of D'MD' that now have positive row sums is at least ^ ^1 — (5 — -^j , 

where for some £ < 1, /3 is the fraction of the diagonal entries of S that are less than £ times the 
average diagonal entry. Indeed we prove in Appendix lAl 

Lemma 4.8. With probability at least |, at most | of the diagonal entries of S are smaller than 
(w) 3 ti mes the average diagonal entry. 

So with probability at least 5 • the fraction of rows with negative row sums in DMD that 
now have positive row sums in D'MD' is at least ^ ^1 — g — | (i§f) 3 ^ > 0. 

Thus, we may conclude that MMatrixSolve is expected to terminate after O (logn) iterations, 
as claimed. □ 



5 Final Remarks 

The reason that our interior-point algorithm currently cannot produce an exact solution to gen- 
eralized flow problems is the dependence of our M-matrix solver on the condition number of the 
matrix, even when approximating in the matrix norm. It would be of interest to eliminate this 
dependence. 

It would also be nice to extend the result to networks with gains. The main obstacle is that the 
resulting linear programs may be ill-conditioned. 
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A Proofs for Section [4] 



Lemma 14. 3L Given vectors vi, . . . ,v n £ IR m and constants a,(3,~f,p E (0, 1), for positive constant 
integer k = kjL(a, (3,~f,p), let R be a k x m matrix with entries chosen independently at random 

from the standard normal distribution, and let W{ = y^Ruj. 
With probability at least p both of the following hold: 



< 0n 



Proof. Let Z,- = ,■> and Z = V] 1 ?-, Zi. 

|| to ,;|| ^i— 1 

Let ri, . . . , Tfc be the rows of R, and let Wij = k^ 1 ^ 2 (rj, v^) be the jth entry of Wi. 

Without loss of generality, we assume that all the V{ are unit vectors. Thus for any given i, 
the expressions k^^wn, . . . , k^^Wik are independent standard normal random variables. So the 
expression 



k k\\ Wi \\ 2 Eti(^i) 



has inverse-chi-square distribution, with mean aim variance jrz^wz^ ■ Therefore, Z has mean 



and variance at most ^^^k-A) ' because 



Var [Z] = Var 



,i=i 



^ Cov [Z^Zj] < a/ Var [Z;] Var [Z,-] = n 2 Var [Z;] = fc 2 n 2 Var 
ij=l ij'=l 



z 

k 



So using Cantelli's inequality, we may conclude that 

Var [Z] 



Pr [Z > (1 + 7 )ra] < 



< 



Var [Z] + (1 + 7 - ~~ 2 + {k — 4)(1 - §) 2 ( 7 - 2 )2 



(1) 



By the same reasoning, has chi-square distribution, with mean k and variance 2k. So using 
Cantelli's inequality, we find that 



Pr [Z < 1 - a] = Pr 



k k 

> 



Zi I- a 



< 



Var [k/Zi 



Var [k/Zt] + (j^-k) 2 2 + {^) 2 k 



Thus, the set {i : Zi < 1 — a} has expected cardinality less than 2 ? 2 nt H. So using Markov's 
inequality, we conclude that 

2 



Pr [|{» : Zi < l-a}\ > f3n] < 



(3(2 + (j^Yk) 



(2) 



Combining inequalities [T] and [5] via the union bound, we find the probability that (i) and (ii) both 
occur is at least 

2 2 

1 



2+ (k — 4)(1 — f)2( 7 - ^)2 0(2 + 
which is greater than p for sufficiently large k. 



□ 
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Lemma 14.21 (Random Scaling Lemma). Given an n x n M-matrix M , and positive real values 
( < 1 and r < j, let D be a random diagonal n x n matrix where each diagonal entry d{ is chosen 
independently and uniformly from the interval (0, 1). 

Let T C [n] be the set of rows of MD with sums at least r times the pre-scaled diagonal, i. e. 

T = {i e [n] : (MDl)i > rm u } 

With probability at least we have 

n*(H)(-'-*)- 

where (3 is the fraction of the diagonal entries of M that are less than £ times the average diagonal 
entry. 

Proof. Let Mo denote the matrix containing only the off-diagonal elements of M. Thus, Mo has 
no positive entries. 

Let B be the set of rows of M in which the diagonal entry is less than ( times the average 
diagonal entry. Thus \B\ = (3n. 

We define a subset J of rows of M whose sums are not too far from being positive. In particular, 
we let J be the set of rows in which the sum of the off-diagonal entries is no less than — § times 
the diagonal entry: 

J=S [ ie[n}: (Mol n )i > —rriir 
Let us prove that J cannot be too small. Let S be the sum of the diagonal entries of M. We have: 
S = l£Ml„- J^(M l n )i 

ie[n] 

> — (Mpl n )j (because M is positive definite) 

ie[n] 

> — ^2 (Mol n )i (because Mo is non-positive) 

ie[n]-(JUB) 

3 
2 



> - ^2 m a (by definition of J) 



iS [n]- (JUS) 

> -— - (J U B)\ (by definition of B) 
2 n 

>~(n-\J\-/3n) 

So we see that \J\ > (1 — (3 — 

Next, let us show that the rows in J have a high probability of being in T. Consider the ith 
row sum of MoD: 

{M Dl)i = ^2d j m ij = - J^my + ^{dj - -)my = -(M l)i + J^(dj - -)m ij 

Since each (dj — g) is symmetrically distributed around zero, we may conclude that |(Mol)i is 
the median value of (MoDl)i. We may also note that (MDl)j = (M(jDl)j + dimn, and that the 
values of (MqDI); and djmjj are independent. 
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We thus have, for i G J: 



Pr [(MDl)i > ruin] > Pr 



{MoDl)i > -(M l)i 



Pr 



2 



> - • Pr 

~ 2 



Pr 



diinu > rmu - -(M 1) 



<U > r + 



d i m ii >rm ii --(M l) 



(by definition of J) 



Thus the expected size of J — T is at most (r + |) \J\. So we find 



Pr 



> Pr 
= Pr 



l^nr|>(i-^)|j| 

Vf-'i"\<[\ + l)\-n 



>i r + * 



r + Z 



(by Markov's inequality) 



1 - 4r 



4r + 7 

The lemma then follows from the lower bound on \J\ proven above. □ 

Lemma 14.61 D' is a positive diagonal matrix. 

Proof. D' 2 = Yi^^Dr is trivially positive diagonal by construction. 
To check that D[ is positive, we use Lemma lA.ll which implies that 

D[l > -M^Mi 2 D' 2 ln-u + 5 (M^l n - V - -X^^ln-v 

To see why the above expression is positive, recall from Fact 14.11 that M-^ 1 and —M\ 2 are 
positive matrices. Furthermore, note that the diagonals of M-Q 

Lemma A.l. 



\D[l - M^{-M X2 D' 2 + < -5\^ ax 



Proof. Recall from the algorithm that 

D^D^l = STSolvepiMn-Di, D X (-M X2 D' 2 + 61)1, e 2 ) 
Therefore, STSolve guarantees that 

\\D^D[1 - D^M U \-M 12 D' 2 + 5I)l\\ DiMiiDi < e 2 \\D^M^\-M 12 D 2 + 5I)l\\ DiMnDi 
or equivalently 



D[l - M^(-M l2 D' 2 + 6I)l\\ Mii < e 2 \\M^(-M 12 D' 2 l + 51) 
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which in turn implies that 

||£>il - M U \-M 12 D' 2 + 6I)l\\ < e 2 K 1/2 \\M^{-M 12 D' 2 + 5I)l\\ 
We can then see 

||.Dil - M^(-M l2 D' 2 + 5I)l\\ < e 2 K X / 2 \\-M^M 12 D' 2 l + tfM^lH 

< e 2 K 1/2 llM^MiaD^lH + 5e 2 K 1/2 \\M^l\\ 

< e 2 K X l 2 \\M^-M 12 D' 2 l\\ + 5e 2 K 1 ' 2 X^ ln n 1 l 2 

< e 2 Kn 1 / 2 \\D' 2 l\\ + Se 2 K 1,2 KknP}-/ 2 (by Lemma EH]) 



< e 2 nn 1/2 

< 2e 2 Kn l/2 



s -i/ 2l 



s D 1/2 i 



+ de 2 K 1/2 X^ n n^ 2 (D' 2 < 5T 1 / 2 by construction) 
+ ( 5e 2 K 1/2 A-J n n 1 / 2 (5T 1 / 2 < 2S~ 1/2 by Lemma M. 



< 2e 2 K\J 2 n + 5e 2 K 1 / 2 X m ] n n 1 / 2 (S D 1/2 1 < \J 2 1 by Fact EM} 



25- 1 e2K 2 A^n + 6 2K 1 /2 n i/2U A 



< -<5A 



-i 



□ 



Lemma 14.81 Wii/i probability at least \, at most \ of the diagonal entries of S are smaller than 



99 \ 3 



.101 



) times the average diagonal entry. 



Proof. Recall that the diagonal entries of S are s~u = where su = II (I — A\ Ai)af\\ 2 and 

a, = \\{R-QD x A x )aff. 

Let us define W{ = jjj \\R(I — M^ 1 A\) af \ \ , where k = kjL (j^, g, ygjj, |)- By Lemma |4T3| 
there is at least | probability that 



1 " o . . 

v ^— ' Wi 

1=1 



< 1.01 



and 



i : — < .99 



1 

< -v 
~ 5 



So, by Lemma [A.2I below, there is at least a | — > | probability that the average diagonal 
entry of S is at most 



1 1 

1 \ -< Sjj _ 1 \ Sjj Wi 
1/ — J rr- 7y ^ — ^ in- rr- 



< 



1.01 



fc(.99) 2 

and similarly we have the following bound on the number of small diagonal entries: 



v £ — ' a. 

i=l 



V Wi Ui 

i=i 



i : — < 

cr,; 



.99 



fe(l.Ol) 2 



1 

< -v 
~ 5 



□ 



Lemma A. 2. With probability at least 1 — it holds for all i that 



Wi 

< — < 



k(1.0l) 2 ~ a ~ A;(.99) 2 
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Proof. We have: 



1/2 _ k l/2 1/2 



= HKfl-QDiAiJof || - \\(R- RAjM^A^af] 
< \\((R - QD^) -(R- RAjM^A^aJW 
= \\(RA^M^A 1 -QD 1 A 1 )af\\ 



N -7=1 



£ Wr-jAjM^D^ - q 3 \\ 2 DlMllDl 

\ 3=1 



\ I' r 3 A l M ll D l 1 ~~ ^'II^MnDi 



< A 1 / 2 

(||aj|| 2 is ith diagonal of M22, so cannot exceed M22's largest eigenvalue) 



A l/2 (_£«U/2 



£ WrjAjM^D^ 1 - qj \\l lMllDl (using FactE© 
\ 3=1 



= («s«) 1/2 
< (^) 1/2 ei 



^Hr^fM^^ 1 \\ 2 DlMllDl (by guarantee of STSolve) 



.005 • s-/ 2 (1.01mn)" 1/2 , 



£||rX^fMl 
\ 3=1 



= .005 • sV 2 (1.01mn)~ 1/2 
■cm-s^k^ilMn)- 1 ' 2 



, \\ r j\\ 2 (because AfMnAi is a projection matrix) 

\ 3=1 



The above inequality does not hold with probability at most based on the fact that expression 
Ej=i Ikill 2 has chi 

-square distribution with mk degrees of freedom. 



< .01 • k 1/2 w] /2 (Lemma [43 implies that su < 1.01 • 
So we conclude that 

< .01 • k 1/2 



nwi 



lf_ _ k l/2 
Wi 



□ 



Lemma 14.71 



MD'l > 



py 2 {SD R -lS D )l 
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Proof. 



MD'l 



MnD[l + M 12 D' 2 l 
M? 2 D[1 + M 22 D' 2 l 

D[l + M^M X2 D' 2 l - SM^l 1 




> 



> 



> 





SD' 2 l 


SD' 2 l 


SD' 2 l 



+ M 



+ 



51 



XmaxWD^l + M u l M 12 D' 2 l - 5M^1\\1 + 
3 



51 

-S\\M^M^1\\1 



:51 + 



51 

-S^nl 



(using Lemmas IA.1I and IA.3D 







> 



> 



SD' 2 l - 25n l / 2 nl 






SE~^D R 1 - \sfSfl 



(using Fact I4.1I5|) 



(using Lemma 



Lemma A. 3. For all positive vectors v , 

\\Ml 2 M^-v\\ < n l/2 n l/2 \\v\\ 
Proof. Define c = A~| n /t _1 / 2 n -1 / 2 \\v\\ = X^ n L ax K 1 / 2 n^ 1 / 2 \\v\\. 

\\Ml 2 M^v\\ < \\Ml 2 M^v\\i = -l T Mj 2 M^v (by Fact ECO M^ 1 and -M 12 are nonnegtr 



= i \ v T M^v + (?1 T M 22 1 - [v T M^ cl T ] M 

< — (v T M^v + c 2 l T M 22 l) 

< \ \ + cXmaxU j = K 1/2 7l 1/2 \\ V \\ 



cl 
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B Solving Matrices from the Interior-Point Method 

In the interior-point algorithm, we need to solve matrices of the form 

M + vv T = 



AD 2 A T + D 2 AD 2 
D\A T D\ + Dl\ 



+ vv 



where A is an n x m matrix with entries bounded by U in absolute value, AA T is an M-matrix, 
and D\,D2,Dz are positive diagonal matrices. We show how to do this using our MMatrixSolve 
algorithm. 

Consider the Schur complement of M: 



M s = (ADfA + D$) - AD((D( + D£)~ D(A = AD(D$(D( + £3) ^ + £> 2 = A S A S 

where As = [ADiD^Df + D 2 )~ l l 2 D2] • Note that Ms is also an M-matrix, and that the eigen- 
values of Ms fall in the range \d 2 min ,d 2 max (JJ ^pmn + 1)] where d m i n and d max are respectively the 
smallest and largest diagonal entry in D\ , D 2 , D3 . 

We can build an solver for systems in M from a solver for systems in Ms, by using the following 
easily verifiable property of the Schur complement: 



Lemma B.l. For M 



M n M12 
Mf 2 M 22 



and Schur complement Ms = Mn — MyiM 22 Mf 2 , we have 





~x{ 


- M- 1 


V 






x 2 









M 



\ Xl - Mg l {bi - M 12 M 2 - 2 1 6 2 )|| Ms + \\x 2 - M 22 \b 2 - M^ Xl )\\ M22 



Then, to solve systems in M + vv , we can use the Sherman-Morrison formula: 

M~ 1 vv T M- 1 



(M + vv T Y l =M _i - 



1 + v T M- 1 v 

In particular, we give the following algorithm, which runs in time O ^mlog ^fcii^ : 



x = Solve(M + vv T , b,e) where M 



AD 2 A T + D\ AD 2 
D\A T D 2 + D' 3 



and b = 


V 


and v = 


Vl 








V2_ 



• Define e 1 = |(1 + v T M- 1 v)- 1 and e 2 = min{±, ^(1 + v T M- l v)- 1 } 

. y> = MMatrixSolve (A s , b x - AD 2 (D 2 + D 2 )~ 1 b 2 , e 1; d 2 min , d 2 max {U^i + 1)) 



y 



• z' 



z = 



{Dl + Dl)-\b 2 -DlA T y') 
MMatrixSolve (A s , «i - AD 2 (D 2 + D 2 )- 1 v 2 , e 2 , d 2 mm , d 2 max (U^ + 1)) 
z' 

(D 2 + D 2 )- l (v2-D 2 A T z')_ 



• Return x = y — _, zzT T b 

a l+v 1 z 
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Lemma B.2. x = Solve(M + vv T ', b, e) satisfii 

\\x — (M + vv T )- 1 b\\ M+vvT < e \\(M + vv T )- 1 b 
Proof. We first show that ||y - M~ l b\\ M < ei j|A^ _1 6 

\\y- M - lb \\ M = 



\M+vv T 



= \\y' -Mg l {b x -AD\{D\ + Dl)- 1 b 2 )\\ Ms (by Lemma EB 
< ei llMg^fti - ADf(Dl + Z?|) _1 6 2 )|| Ms (guaranteed by MMa 



^ifll^&IL- 



|M 22 1 6 2 | 



m 22 



< ei llAf" 1 ^! 



guaranteed by MMatrixSolve) 
(by Lemma iB.lj) 



ill 



By the same reasoning, ||z — M 1 ^|| A// < £2 -M a^- 

Next, let us define the inner product (v\, v 2 ) M = v\Mv 2 . We will use repeatedly the inequality 
\(vi,v 2 ) M \ < \\vi\\ M \\v 2 \\ M 

Recall that we return the value x = y — ** v t z ■ So we begin by analyzing the expressions zz T b 
and v T z: 



\zz T b - M~ 1 vv T M~ 1 b\\ M < 



< 



< 



zz T b - zv T M- 1 b\\ M + \\zv T M~ 1 b - M- l vv T M- 1 b\\ M 
{z-M- l v,M~ x b) M \ Hzll^ + |<M- 1 »,M- 1 6) M | \\z-M~ 1 v\ 



M 



z ~ m ~ 1v \\m \\ M ~ lb \\M Mm + \\ m ~ 1v \\m \\ M ~ lb \\M \\ z 



z-M^v\\ M \\M-H\\ M 



zl^ + IlM- 1 ^! 



M) 



z - M^v\\ M WM-H^ (\\z - M-\\\ M + 2 llM-^HJ 



<e 2 (e 2 + 2)||M- 1 t;|| 2 M ||M- 1 6 



M~ 



(3) 



\v T z - v T M~ 1 v\ = 



v\ = \{M~ l 
< WM^v 



v,z -M" 1 



v 



M\ 

m\\ z ~ M ~ 1v \ 
= e 2 (v T M~ l v) 



M 



(4) 



We thus have: 
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\x- (M + w T ) _1 6| 
zz T b 



M 



(5) 



1 + V 1 z 



< 



< 



< 



< 



y-M~ 1 b\\ M + 



, M- 1 vv T M- l b 

M~ l b = , 

1 + v T M _1 v 

zz T b M- 1 vv T M- 1 b 



M 



1 + V T Z 

1 

• !/ 1 + V T Z 

1 



1 + V T Z 



+ 



M 



M~ 1 vv T M- 1 b M- 1 vv T M- 1 b 



1 + v T z 



1 + v T M~ l 



v 



M 



\zz T b - Ar 1 vv T M- 1 b\\ nr + \ vTz ~f M lv \ \\ M - 1 vv T M- 1 b\ 
1 llM l + v T M~ 1 v 11 1 



M V T Z 



y- M ~ lb \\ M + 
y- M ~ lb \L + 

y- M_lb llM+ (i- e2)v T M -x v 
(by equations [3] and E|) 

e 2 (e 2 + 2)||M- 1 - i ;||' / ||M- 1 6|| M + e 2 |<M- 1 t;,M- 1 6> M | \\M~ 1 v\ 



\zz T b - M- l vv T M- l b\\ M + \— rp ^ M ^ \\M- 1 vv T M- 1 b\\ M 
1 



— (e 2 (e 2 + 2) llM-^H^ \\M~ 1 b\\ M + e 2 HAr 1 ™ 2 ^-^ 



Ai 



y- M ~ lb \\ M + 



M 



y-M-'b\\ M + 



(i-6 2 )HM-^r M 

, ;J| , 6 2 (e 2 + 2) \\M^vf M \\M^b\\ M + e 2 ||M-^||^ \\M^b\\ M 



(l-e 2 )\\M- 1 v 



1 l|2 



1 — e 2 



< 



lAf 

( £1 + ^)||M-.,| 
< e(l + u^- 1 ^)- 1 ||Af -1 &| 
So we conclude 



M 



\x — (M + ^■y T )~ 1 b|| M+ ^ T < (1 + ^ T M-^) 1/2 ||z - (M + ui; T )" 1 6 

< e(l + v T M- l V y 1 / 2 ||M _1 6|| M 
= e{l + v T M- l v)- 1 ' 2 \\b\\ M . 1 
^ 6 ll & ll(A//+™r)-i by LemmaEStii 
= e ||(M + ^)- 1 6|| MW 



Lemma B.3. For all vectors v , w, and symmetric positive definite M : 



(6) 



M by Lemma fB.3n ) 
by equation ^ 



□ 



(i) 
(n) 



\w\ 



M+vv 1 



W 



(M+vv T Y 



< \\w\\ M (1 + v T M~ x v) 1/2 
x > Hwll^-iCl + ^Af- 1 !;)- 1 / 2 
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Proof of (i). 



Hm+wT = (w T (M + vv T )w) 1 / 2 
= (w T Mw + {w T vf) l/2 
= (\\w\\ 2 M + (w,M- 1 v) 2 M )^ 
^(IkllM + ll^llLllM-^l^) 1 / 2 
= \\w\\ M {l + \\M^vf M f/ 2 



Proof of (ii). 



w 



\ {m+vv t ) - 1 = (w t (M + vv t )- 1 w) 1 / 2 



> 



w T M" 1 



w T M^ 1 w 



M- 1 vv T M~ 1 \ \ 1/2 
w 



1 + /M- 1 )) J J 
w T M~ 1 vv T M~ 1 w 



1 + v T M~ 1 v 



1/2 



= \\w\ 



(w,v) 



1/2 



\w\ 



M- 



M- 



M- 



1 + \\v\ 



M- 1 



\W\ 



M 



-i \\V\ 



1/2 



M- 



i + IM 



M~ 



1/2 



VI 



\w\ 



M- 1 



M- 



i + IMlW 

1/2 



\w\\ M -i (l + \\v\\ 2 M -i) 



27 



C Interior-Point Method using an Approximate Solver 

Throughout this section, we take Solve to be an algorithm such that x = Solve(M, b, e) satisfies 

||x-M -1 6|| M < e||M -1 6|| 

We use the notational convention that S denotes the diagonal matrix whose diagonal is s. The 
same applies for X and x, etc. 

lfc denotes the all-ones vector of length k. 

(l denotes the interior of polytope f2. 

We are given a canonical primal linear program 

z* = min { c T x : Ax = b; x > 0) 

X 

which has the same solution as the dual linear program 

z* = max { b T y : A T y + s = c; s > 0} 

where iisannxm matrix, x, s, c are length m, and y, b are length n, and m > n. (Unfortunately, 
this use of n and m is reversed from the standard linear programming convention. We do this to 
be consistent with the standard graph-theory convention that we use throughout the paper.) 
We let Q D denote the dual polytope 

n D = {{y, s) : A T y + s = c; s > 0} 

so we can write the solution to the linear program as z* = maXf y s \^QD b T y. 

In this appendix, we present an InteriorPoint algorithm based on that of Renegar [Ren88], 
modified to use an approximate solver. Our analysis follows that found in [Ye97]. 

Theorem 12.11 x = InteriorPoint(^4, b, c, \ m in, T, y , e) takes input that satisfy 

• A is an n x m matrix; b is a length n vector; c is a length m vector 

• AA T is positive definite, and \ m in > is a lower bound on the eigenvalues of AA T 

• T > is an upper bound on the absolute values of the dual coordinates, i.e. 

II f Hoo < an d \\ s \\od < T f or a M {v> s ) th a t sa ti s fy s = c — A T y > 

• initial point y° is a length n vector where A T y° < c 

• error parameter e satisfies < e < 1 

and returns x > satisfying \\Ax — b\\ < e and z* < c T x < z* + e. 
Let us define 

• U is the largest absolute value of any entry in A,b, c 

• s^j n is the smallest entry of s° = c — A T y° 

Then the algorithm makes O ( ^/m log TU ™ — ) calls to the approximate solver, of the form 

Solve (AS~ 2 A T + vv T , -, e') 
where S is a positive diagonal matrix with condition number O ( T ^ m ^ ; and v,e' satisfy 

, II ^ II m(y TUm \ 
log — = (log -o—J 
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C.l The Analytic Center 

Standard interior-point methods focus on a particular point in the interior of the dual polytope. 
This point, called the analytic center, is the point that maximizes the product of the slacks, i.e. 
the product of the elements of s. For the purpose of our analysis, we use the following equivalent 
definition of the analytic center: 

Fact C.l (see |Ye97l §3.1]). The analytic center of Q D = {(y, s) : A T y + s = c; s > 0} is the 

unique point (y, s) 6 Cl D that satisfies rj A (s) = 0, where we define 

x A (s) = ST 1 {I - S- 1 A T {AS- 2 A T )~ 1 AS~ 1 )l m 

VA (s) = \\Sx A (s) - l m || = ||S'- 1 ^ T (A5- 2 A r )- 1 A < S- 1 l m || = WAS- 1 !^^^ 
These definitions of x A and r\ A satisfying the following properties: 
Lemma C.2. Let (y, s) be the analytic center of Q. D . For any point (y, s) G Cl D we have 

(i) Ax A (s) = 

(ii) r] A (s) < 1 implies x A (s) > 
(Hi) x A (s) = S~ l l m 

(iv) For all x satisfying Ax = 0, it holds that \\Sx — l m || > r] A (s) 

The first three properties are straightforward from the definition. We present a proof of the 
last: 

Proof of \ C. glf ivli Note that Sx A (s) — l m is orthogonal to S(x — x A (s)), because 

(Sx A (s) - l m ,S(x - x A (s))) = (S- 1 A T (AS- 2 A T )- 1 AS- 1 l m ,S(x - x A (s))) 

= ((AS- 2 A T )- 1 AST 1 l TO ,A(a; - x A (s))) 
= <(A5- 2 A T )- 1 A5- 1 l m ,0> 
= 

We thus have 

\\Sx - l m \\ = \\Sx A (s) -l m + S(x - x A (s))\\ > \\Sx A (s) - l m \\ = r] A (s) 

□ 

It will be useful to note that the slacks of the analytic center cannot be too small. We can 
bound the slacks of the analytic center away from zero as follows: 

Lemma C.3 (compare |Ye97t Thm 2.6]). Let (y,s) be the analytic center of ft D . For every 
(y, s) G VL D , we have s > 

Proof. 



S 



-l . 



= m + l'^ n S' l {s - s) 

= m + llS- 1 {(c-A T y)-(c-A T y)) 

= m + l T J- 1 A T (y-y) 

= m 

where we know from Lemmas IC.2fi)l and IC.2nii)1 that AS~ l l m = □ 
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y+ = NewtonStep(A, c, y 

• Let s = c — A T y 

• Let d y = Solve (AS* -2 ^, -AS~ l l m , e 3 ) where e 3 



Return y + = y + (1 - e 3 )d 2/ 



Figure 2: Procedure for stepping closer to the analytic center 

Let us note that a point (y, s) G £l D that satisfies rjA(s) < 1 is close to the analytic center, in 
the sense that the slacks s are bounded by a constant ratio from the slacks of the analytic center: 

Lemma C.4 ( |Ye97[ Thm 3.2(iv)]). Suppose (y, s) G il D satisfies t]a(s) = T) < 1 and let (y, s) be 
the analytic center of Q D . Then — l m || < j^j- 

If (y, s) G (l D is sufficiently close to the analytic center (as measured by tja), then with a single 
call to the approximate solver, we can take a Newton-type step to find a point even closer to the 
analytic center. This MewtonStep procedure is presented in Figure [2 

In the first part of the following lemma, we prove that the point returned by NewtonStep is 
indeed still inside the dual polytope. In the second part, we show how close the new point is to the 
analytic center: 

Lemma C.5 (compare [Ye97[ Thm 3.3]). Suppose (y, s) G tl D satisfies t]a(s) = tj <1. 
Let y + = NewtonStep(A c, y) and s + = c — A T y + 

Then (i) s + > and (ii) tja(s + ) < rj 2 + ^rj 

Proof, (i) The solver guarantees that 

\\d y + (^ 2 ^ T )~ 1 ^ 1 l|| yl5 _ 2AT < e 3 || (AS- 2 ^)- 1 ^- 1 !!!^.^ = e 3 • r, 
or equivalently 

\\S- 1 A T d y + S- 1 A T (AS- 2 A T y 1 AS- 1 l\\ < e 3 ||S _1 A r (AS ,_2 A r )- 1 AS _1 l|| = e 3 • rj (7) 

and so 

\\S~ 1 A T d y \\ < (1 + e 3 ) \\S- 1 A T (AS- 2 A T )- 1 AS- 1 1\\ = (1 + e 3 ) ?? < 1 + e 3 
We thus have 

\\S- 1 s + - l|| = 1 1 ^ 1 (s - (1 - e 3 )A T d y ) - l|| 
= (l-e 3 )\\S- 1 A T d y \\ 

< (l-63)(l + 63) < 1 

Thus S~ x s + is positive and so is s + . 
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(ii) Let x = xa{s) and x + = xa{s + ). We have 

Va(s + ) < \\Xs + - l m || (by Lemma [021w)J) 
= \\X(s-(l-e 3 )A T d y )-l m \\ 

= ||(1 - e 3 )IS(I S - l m - S~ 1 A T d y ) - (1 - e 3 )(XS - - l m ) + e 3 (I S - l m )|| 

< (1 - e 3 ) ||X5(5- 1 A T d 1/ - Sx + l m )|| + (1 - e 3 ) \\{XS - I)(Xs - l m )|| + e 3 - l m )|| 

< (1 - e 3 ) ||5aj|| ||S- 1 A T d 2/ - Sx + l m || + (1 - e 3 ) ||5aj - l m || 2 + e 3 ||5a; - l m || 

(using the relation ||Vtw|| < \\v \\ \\w\\ < \\v\\ \\w\\) 

< (1 - e 3 )(||5a! - l m || + ||l m ||) \\S- 1 A T d v -Sx + l m \\ + (1 - e 3 ) \\Sx - l m \\ 2 + e 3 \\Sx - 1 
= (1 - e 3 )(r7 + y/m) \\S- l A T d y - Sx + l m || + (1 - e 3 )7? 2 + e 3 r/ 

= (1 - e3 )(r? + y/m) \\S- 1 A T d y + 5- 1 ^ T (A5- 2 ^ T )- 1 A^ 1 l m || + (1 - e 3 )r/ 2 + 

< (1 — e 3 )(ry + \/m)e s rj + (1 - e 3 )?? 2 + e 3 r/ (by equation [7]) 

< e 3 (ry + v 7 ™)?? + (1 - + 
= rj 2 + e 3 (y / m + l)r? 



C.2 The Path-Following Algorithm 

In a path-following algorithm, we modify the dual polytope Q D = {(y, s) : A T y + s = c; s > 0} 
by adding an additional contraint b T y > z, where z < z*. As we let z approach z* , the center of 
the polytope approaches the solution to the dual linear program. 



Letting s 



gap 



b y — z denote the new slack variable, we define the modified polytope: 



(y,s,s 



gap) 



A T y + s 
b T y + Sg a p 



Sgap ^ 



Using a trick of Renegar, when we define the analytic center of , we consider there to be m 
copies of the slack s gap , as follows: 



Definition C.6. The analytic center of is the point (y,s,s gap ) £ Cl^ ' z , that satisfies 

fj(s,s gap ) = ; where we define 

V(s,Sgap) = ri A (s) 



where A = [A —61^] and s 



^gap^-m 



The central path is the set of analytic centers of the polytopes < ^l^ z \ 

A path-following algorithm steps through a sequence of points near the central path, as z 
increases towards z* . It is useful to note that given any point on the central path, we may easily 
construct a feasible primal solution follows: 



gap) 



be the analytic center of £lP z - 



Then the vector x = ^5 _1 l r 



Lemma C.7. Let (y, s, s 

satisfies Ax = b. More generally, for any (y, s,s gap ) E ^ z , the vector x = ^ £ 5~ 1 l m satisfies 



\Ax 



m 



V{ s i s gap) 
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s+, sj ap , z+) = Shif t(y, s, s gap , z) 



Let z+ = z + ^E- 



Let t/ + = NewtonStep(A, c, y) where = [A —61^] and c 



c 

-z+l r 



Let 







' c - A T y + ' 


L gap. 




b T y+ - z + _ 



Figure 3: Procedure for taking a step along the central path 



Proof. We prove the second assertion: 

\ r - h\\ -- - - * 9ap II /W 1 ! -mr'/il 
• lx °H(A5- 2 A T )- 1 — m ll^ 10 J " m mS Sap°ll(A§- 2 A T )- 1 

(AS~ 2 A T ) 



_ s gap 1 


\AS' 


l l - 


m ' 




Sgap 


AS' 


J- 2m 


m 






_ s gap 

m 


rj(s, 


Sgap) 



The first assertion now follows from the definition of analytic center. □ 

Let us now describe how to take steps along the central path using our approximate solver. 
In Figure El we present the procedure Shift, which takes as input a value z < z* and a point 
(y, s, Sg a p) G Cl® z satisfying r](s, s gap ) < j^. The output is a new value z + that is closer to z*, and 

a new point (y + , s + , s^ ap ) G &® z + satisfying rj(s + , s^ ap ) < jq. The procedure requires a single call 
to the solver. 

Let us examine this procedure more closely. After defining the incremented value z + , if we 
let s' gap = b T y — z + = s gap — (z + — z), then (y,s,s' gap ) is a point in the shifted polytope ^ z +- 
However this point may be slightly farther away from the central path. One call to the NewtonStep 
procedure suffices to obtain a new point (y + ,s + ,s+ ap ) G fl® ' z+ that is sufficiently close to the 
central path, satisfying r/(s + , s+ ap ) < jq. 

We prove this formally: 

Lemma C. 8 (compare |Ye97l Lem 4.5]). Given z < z* and (y, s,s gap ) G Cl® z where fj(s,s gap ) < j^, 
let s 'gap = bT y ~ z+ an d s+ i sf ap , z + ) = Shift(y, s, s gap , z). Then 

(i) z + < z* 

(ii) s' g a P > and fj(s,s' gap ) < ^ 
(Hi) s + ,sj ap > and f]{s + ,sj ap ) < i 
Proof, (i) 

z+ = z + b ^~ Z <z + {b T y -z) = b T y < z* 
lU\/m 



(ii) We note that s' gap = s gap - (z + - z) = (l - j^^j s gap > 0. 
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Let us write s 



^gap^-m 



and s' 



s' 1 



(S 



gap s gap)^-r 



and note that 


Z + -Z)l r 










Sgap -f 




lOv/m m 



(8) 



Let us define 5 



x 



x 



yap. 



xx(s). So we have 



fj(s,s' 



gap) 



< 
< 



< 
< 



S'x - 


" l2m 


Sx - 


l2m 


Sx - 


l2m 


Sx - 


l2m 


Sx - 


l2m 


Sx - 


l2m 


Sx - 


l2m 


11 ~( 
10 A ' 


Sgap) 



(by Lemma OExB 

(5 -5)5 
1 



ioV^ l|S9ap:C9apl1 

1 



(by Equation [8|) 



1 



ioV^ l|S9ap:C9ap im " ^ 10 

1 II - II j 1 



Imll + 



1 



+ 



1 

To 
l 

To 

11 11 21 

- To ' To + To ~ Too 



55-1 



2m 



+ 



1 

To 



0) 



(iii) By Lemma fC. 51 we have s + , s+ ap > and 

. i I . , . , 2 1 , , , / 21 \ 2 1 21 1 
•7(- ,4a P )<V(s,s gap ) +^v(s,s 9ap )<{- m ) + 20'T00 < T0 

□ 



We now present the complete path-following InteriorPoint algorithm, implemented using an 
approximate solver, in Figured! For now we postpone describing the FindCentralPath subroutine, 
which gives an initial point near the central path. In particular, it produces a z c < z* and 
(y c ,8 ,Sg ap ) E ) c satisfying fj(s , Sg ap ) < j^. Once we have this initial central path point, 
Lemma IC.81 tells us that after each call to Shift we have a new value z < z* and new central path 
point (y,s,s gap ) £ £1® Z that satisfies fj(s,s gap ) < jq. 

Later we will analyze the number of calls to Shift before the algortihm terminates. First, let 
us confirm that the algorithm returns the correct output: 

Lemma C.9. The output of x = InteriorPoint(A, &, c, y°,e) satisfies 

(i) x > 

(ii) \\Ax — b\\ < — t=£- -. ,„ 

1/11 11 — 12 v / 2-Tn 1 /2 

(Hi) c T x < z* + e 
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x = InteriorPoint(^4, 6, c, y°, e) 

• Compute (y c , z c ) = FindCentralPath(^4, 6, c, y°) and 

• Set (y, s, s gap , z) := (y c , s c , s^ apl z c ) 

• While Sgap > §: 

- Set (y, s, s ga p, z) := Shif t(y, s, s 3ap , z) 

• Compute t> = Solve(AS~ 2 A T , ASHl 2m ,e 4 ) 

where 1 = [A -61^] 
and s 



c - A T y c 
& T y c - z c 



and e 4 = min 1 



• Return x 



Sgap^-m_ 



where 



TU n I 

and s m in is the smallest entry of s 



' x' ' 




x' 

L Uj gap_ 





S- 1 ^ - S~ 2 A T V 
b gap * a gap u v . 



Figure 4: Dual path-following interior-point algorithm using an approximate solver 



Proof, (i) To assist in our proof, let us define x' 
We have 



x 

L x 'gaplm 



and note that x' = S 1 l 2m -S' 2 A T v. 



Sx'-l 



2 m 



S- l A T v 
v \\as- 2 a t 



< 



{AS- 2 A T )- 1 A~S- 1 l 2m 



AS- 2 A T 
-2 aT\-1 15-1 



+ 



v - {AS^A^AS' 1 ! 



2m 



< (l + e 4 ) (AS-'A'y'AS-'l 



2m 



AS' 2 A T 



AS- 2 A T 
(by guarantee of Solve) 



= (1 + e 4 ) • ?y(s, Sgap) < 2 • s ga p) < 2 • — < - 

Since s is positive, we conclude that x' must also be positive, and so must be x. 
(ii) We have 

1 



(10) 



\Ax - 61 



mx 



, u Ax' - mx gap b\ 

gap 



1 



mx 



gap 



mx 



gap 



Ax' 

AS- l l 2m - AS~ 2 A T v 
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Observe that the largest eigenvalue of the matrix AA T = AA T + mbb T is less than the trace, which 
is at most 2nmU 2 . Thus, the largest eigenvalue of AS~ 2 A T is at most 2nmU 2 s^ r 2 n . So we proceed 



< 



{2nmU 2 s- 2 y/ 2 



< 



< 



< 



mx gap 

{2n) 1 l 2 U 

,1/2™/ 

■ L gap- 
(2n) 1 / 2 U 

,1/2™/ 

^gap"- 

{2n) l l 2 U 

Tfl^/ 2 Xg a pS m i n 

{2n) 1 / 2 U 

,1/2, 

•^gap"- 

(2n) l l 2 U 

Vn}'l 2 X gapSmin 
1 



AS~ 2 A T v 



AS 1 l2m 

(^- 2 i T )- 1 i5~ 1 l 2m - < 
e 4 |(i5- 2 i T )- 1 i5- 1 l 2m 

e 4 • fj(s,Sgap) 



(AS- 2 A T )- 



AS- 2 A T 



AS~ 2 A T 



1/2 

nTU 



1 

10 



1 



5^ • TnV2 x ' 



< 
< 
< 



1 



!!"!> 



■ S, 



(by guarantee of Solve) 



5^2 -TnV2 4 

1 5 e 

5^ • TnV2 ' 4 ' 3 
e 



(from equation 1 101 we know s 



:iap x gap — 



12^/2 -Tn 1 / 2 



(11) 



(iii) We have 



T 
S X 



T 1 
S X 



mx gap 



5 x^ 6 
> — • Sg a p (from equation [101 we know s ga px' gap < -) 

0772 O 

5115a'!!! 

6m gav 

5 

— (Hl m Hl ~~ ~ ' S 9 a P 



6m 
5 



(m — lliS'x' — lm||i) 



6m 
5 



1 J ' s gap 



> — (m - VmllSV - l m ||) 
6m 



s 



gap 



> ^m - ^V^^ • s gap (by equation [TP]) 



5 4 

* * 777 * S oclu 

6m 5 



'gap 



(12) 
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T 
S X 



T l 
8 X 



mx 'gap 



5s T x' 4 
< ^ ■ s gap (from equation [101 we know s gap x' gap > -) 

_ 5jSx% 
~ 4m ' Sgap 

— [\\^ X ~ -'-"illl + l|l"i|li) • Sg a p 



5 



— (\\Sx' -lml^ + m) 



'gap 



5 

< — (y/m\\Sx' -l m \\+m) ■ Sg ap 

< (- \pm + m^j ■ Sgap (by equation [TO 

5 6 
- 4m" ' 5 ' m ' Sgav 



*gap (13) 



We then have 

c T x — z* > c T x — b T y 

(c T - y T A)x + y T (Ax - b) 



2 

= s T x + y T (Ax - b) > - ■ s gap + y T {Ax - b) (by Equation [T3]) 

3 

< ■e + y T (Ax - b) 
<l. e +\\y\\\\Ax-b\\ 

< - -e + Tn 1/2 \\Ax - b\\ 
3 

2 l 

< — • e H — • e (by Lemma l"C.9f ii)) 

3 12 v 2 

< e 
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T * T 

c x — Z < C X — Z 



(c 1 - y 1 A)x + y 1 {Ax - b) + b 1 y 

s gap 



s T x + y T (Ax - b) + s 6 



< - • Sgap + y T (Ax - b) (by Equation H3 

< --e + y T (Ax - b) 

5 , „ „ 

< 7 • e + y \\Ax - b\\ 
o 

< - -e + Tn 1 / 2 \\Ax - b\\ 

5 \ 

< - • e H -= ■ e (by Lemma fCT9l fii) ) 

6 12\/2 

< e 



□ 



Next, we analyze the number of Shift iterations until the algorithm terminates. We can 
measure the progress of the algorithm with the potential function B{z): 

in 

B(z) = 1°§ *j + ^log Sg a p where (y, s, s gap ) is the analytic center of £1® Z 

3=1 

Soon, we will show how a decrease in B(z) implies that s gap is decreasing and thus the algorithm 
is making progress. Let us first show that the value of B{z) decreases by £l(^/m) after each iteration. 

Lemma C.10 (compare |Ye971 Lem 4.6]). Given (y,s,s gap ) E Cl® z satisfiying fi{s,s gap ) < jq, let 
(y+ z + ) = Shift(y, z). Then B(z + ) < B(z) - @(^/m). 

Proof. Let (y, s, s gap ) and (y + , s + , s^ ap ) respectively be the analytic centers of ttj? z an< ^ z+ • 
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Following Lemma lUTf! we define x = ^f-S~ l l m that satisfies Ax = 6. We have 



B(z + )-B{z) 

e 2m 



u' =1 



'gap 



'gap 



< 



/m ^ — ^ fl ■ 



.gap 



, I J_ \ ^ °J_ + 1 •>'/< ° .'/«/• 



1*2 



'gap 



1 + 



1 



2 J 



gap 
1 



(0 



s) + (s^" ap S 5aj , 



)) 



1 + ( ((c - ^ T £ + ) - (c - A T y)) T a: + (s gQ;) .„„„ 



'gap 
1 



1 + iy^— ((v - & + ) T ^z + ( 



'gap 
1 



s gap s gap 



)) 



1 + ^pr— ((if - £ + fb + (S+p - Igap)) (by Lemma E2D 

1 -.T. 



1 + 



2s fi 



((& T y 



'gap/ 



'gap 

1 + — (z - Z + ) 



gap 



)) 



2s 



<i»p 



?gap 



■ Sg a p 



< 1 



9 



200V^ 

9 



( Sgap < 10 ^ LemmaE ^ 



'gap 



< g 200^m 

We conclude that B(z + ) — B(z) < —^my/m- 



(14) 



□ 



Let us now show that a decrease in the potential function B{z) implies a decrease in the value 



of s 



gap- 



Lemma C.ll (compare |Ye971 Prop 4.2]). Given (y,s,s gap ) G Clf z and (y + ,s + ,s+ ap ) G ' z+ 
where z + > z and fj(s, s gap ) < rj and fj(s + , ) < rj for r\ < 1. Then 



igap 



'gap 



<(l-2r,)-(e 



B(z 1 )-B(z 2 ) 



Proof. Let (y, s, s gap ) and (y + , s + , s+ ap ) respectively be the analytic centers of Q® z and £lf z +- We 
define x = ^S^lm and s + = ^(S" + )- 1 l m , which by Lemma E3 satisfy As = 6 = 
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We have 




So 



l + ^((c-A T y) AT!.+"*~+\ 

Sgap / °gap 

l + ^{y-y + ) T Ax 

jap 



°gap / °gap 

v ■sgap / ^9 a p 



fl + 4- ((c - - (c - A T $+))' J x 1 ■"■*"" 

\ ■Spap / ■'sip 

l + ^(S-S + ) T aj > ) • |^ 

■Sgap / "Sgap 

< \ 1+ _L* a T x \ *9ap 

1 + 



gap / ogap 
gap \ Sgap 



<H + i 



gap / a gap 

, 2 
gap ' 



'gap 



'gap 



B(z)-B(z+) 

> e 2m 



6gap 

Using Lemma IC.41 we may conclude 



* *+ 1 — 71 / 

Sgap ^gap 'Sgap ^gap ^ 1 — T7 / B(z)-B(z ) 



Sgap ^<?ap ^9 a P Sgap 



1— jj 



> -tP- ■ e 2^ - 1 



□ 



Corollary C.12. The InteriorPoint algorithm makes O (y^log^^) calls to Shift 
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Proof. Recall that the algorithm will terminate only when the value of s gap has decreased from its 
initial value of s^ ap to below |. Thus, Lemma fC 111 ensures us that s gap will be smaller than | once 

B(z) has decreased by SI (^m log . According to Lemma fC.lOl this occurs after O (^Jm log ^m^J 
Shift iterations. □ 



C.3 Finding the Central Path 

It remains for us to describe how to initialize the path-following algorithm by finding a point 
near the central path. Essentially, this is accomplished by running the path-following algorithm 
in reverse. Instead of stepping towards the optimum given by b, we step away from the optimum 
given by the vector b = A(S°)~ 1 l m that depends on our initial feasible point (y°, s°) G fl D . 

Our analysis parallels that in the previous section. The following function rj measures the 
proximity of a point (y, s,s gap ) G Cl^ z to the central path given by b: 



v( s ,§ga P ) = 77a(3) where A = [A -61^] and s 



?gaplm 



To initialize the algorithm, we observe that (y°, s°,m) G ^f 2 o is on the b central path, where we 



define z = b y 



T.,0 



m: 



Lemma C.13. r](s°, rn) = 



Proof. Defining s 

Ms )- 1 ! 



s 
ml 



o 



we have 



2;/f 



[A 



m~ l l m 



A(S' 



0\-l- 



1/1 

m 



Thus, rj(s°,m) = AiS )' 1 ! 



*-2m 



{A(s°)- 2 A T y 



□ 



We present the FindCentralPath algorithm in Figure [5j Starting with z = z°, we take steps 
along the b central path, decreasing z until it is sufficiently small that the analytic center of 

is close to the analytic center of O , and therfore also close to the analytic center of ' for some 
sufficiently small z. 

Let us show that the Unshif t procedure indeed takes steps near the b central path: 
Lemma C.14 (compare Lemmas IC.8p . Given (y,s,s gap ) G Q,® z satisfying rj(s,s gap ) < ^j. Let 



-gap 



b y — z + and (y + , s + ,s 



9W 



Unshif t(y, s, s gap , z) . Then 



< 



51 
400 



(*) V(s,2gap) 

(ii) rj(s+,s+ ap ) < ^. 

Proof of \0.1JJ( i) . Following the proof of Lemma l(J.8f i) through equation [9j we have 
U( s ^'ga P ) 



11 , 1 11 1 1 

< ■ 77(S, S oaD ) H < • i 



51 
400 



□ 



40 



(y,s 


Sgapi 


FindCentralPath(,4, b, c, y ^ 








Define 


















" s° " 

s° 




"c - A T y°~ 

m 






i T n 


• 


Set (y, s,s gap ,z) := 


(y°, s°, s^ ap , z ) 






• 


While s 9ap <40Aj/ 2 Tm||6||: 










Set (y, s, Sgapi 


z) := Unshift (y, s, 


—gap i — ) 




• 


Return (y 
and z = b 




40A 


~ 1 ^ 2 Tm \\b 

mm x '"■ II u 










and s gci p = 


= b T y - z 









(y + ,s + ,s+p,£ + ) = Unshift(y,s,s 9ap ,2;) 

• Let z+ = z - -|^= 

lOvm 

• Let y + = NewtonStep(A, e, y) where 4 = [^4 
c - ^y" 4 



-61^] and c 



• Let 



-gap. 



b T y+ - ^ 



Figure 5: Algorithm for finding point near central path given feasible interior point 
Proof of \C.14\ ii)- By Lemma |C. 51 we have 

V(s + ,s_; ap ) < fj(s,s' gap f + ±fj(s,s> gap ) < (JL^ + 1 • |L < _L 

□ 

Next, let us prove that the point returned by FindCentralPath is indeed near the original 
central path (i.e. the path given by &): 

Lemma C.15. For y° satisfying A T y° < c, let (y, s, s gap , z) = FindCentralPath(A, b, c, y°). 

Then (y, s, s gap ) £ ttf z and fj(s, s gap ) < jq. 



Proof. Using the values at the end of the algorithm, we write s 



s 


and 3 = 


s 






Sgaplm_ 




Sgap 1 m 
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To begin, we note 



s gap >m\J*Tm\\b\ 



= mm{T- 2 \ min )- 1 / 2 \\b\ 
> 40m||6|| (As ._ 2A T)-i 



(15) 



where the last inequality follows because the smallest eigenvalue of AS 2 A T is at least T 2 A, 
Similarly, 



40\J/ n 2 Tm\\b\\ > 40m || 6 1 



{as- 2 a t y 



(16) 



We have 

fj{ s i s gap) 



AS-H 2 m 



(AS~ 2 A T )~ 



< 



AS' 1 ! 



2m 



(AS~ 2 A T )~ 1 

(because AS~ 2 A T — AS~ 2 A T = ms~ 2 p bb T is positive semidefinite) 



AS 1 1-2m-ms g a p b + msJ„b 



-l 

3 gap 1 



gap 

(A5- 2 A T )-! 



< 
< 



AS' 1 ! 



AS' 1 ! 



2m 
2m 



(AS~ 2 A T )~ 



+ mSgap || 6 1| (AS- 2 AT)- 1 + ms gap II b II (AS' 2 AT)- 1 



+ h — (by equations [15] and [TBI) 

(AS^AT)- 1 40 40 v — — v 



< ( 1 + ms gap || b \\(AS- 2 A T )- 1 



1/2 



2 m 



1 1 

(il" 2 ^)- 1 40 40 

2 IT /i o-2 4 T -2 i> T\ 



(by Lemma El using the fact that AS~ 2 A T - AS~ 2 A T = ms~ 2 p bb T ) 



1 ' ,'„ 



1/2 






4(5 1 l2m 



1 1 I— 

+ - h — (by equation 15 

(A5-2A^)-i 40 40 v J H 



40 



i+— J -??(s,s 9ap ) + — + 



1 1 
4~0 + 40 



1 

20 



< 2 • rj(s,s ga p) + 

< 2 • ^ + ^ (by Lemma EHii) ) 

1 

~ 10 



To measure the progress of the FindCentralPath algorithm, we define B(z): 

rn 

where (fir, s,s gap ) is the analytic center of Q.® z 



B(z) = ^2 l°g *j + "mlogs 

3=1 



gap 



□ 



Lemma C.16 (compare Lemma [C.lOj) . Given (y,s,s gap ) G Q,^ z satisfying rj(s,s ga p) < ^. Let 
(y+,s+,s+ p ,z+) =Unshift(y, s,s 9ap ,z). T/ien B{z+) > B{z) + Q(<fm). 
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Proof. We will follow the proof of Lemma |C. 101 with some minor changes. 
Before we proceed, let us recall the definition z + = z — to note that 

§'gap = Sgap + Z — Z + = W\/m(z — Z + ) + Z — Z + < lly/m(z — Z + ) (17) 

Now, we switch the places of z and z + , and follow the proof of Lemma lC.lOl up to Equation 1141 

B(z)-B(z+) 1 
e 2m < 1-| — _ . ( Z + - z) 

^Sgap 



We continue: 



s' 

< 1 ^ «+ ( b y Equation [IT]) 

22-y/m • Sgap 

1 349 

< 1 = (by Lemmas ICTTl and [Oil) 

22ym 400 

169 



< g 4400-vAS 



□ 



Corollary C.17. T/ie FindCentralPath algorithm makes O (y/m log . "^T^ — ) caZZs to Unshift, 
where s^m ^ s ^ e smallest entry of s° = c — A T y°. 

Proof. Recall that the algorithm will terminate only when the value of s gap has increased from its 
initial value of m to at least 40A mi ^ Tm \\ b ||. So, by Lemma IC.1H this will have happened once 
B{z) has increased by ^mlog ^A m ^ 2 T ||6||^ J . 

According to Lemma |C.16| this occurs after O (^/m\og ^A m ^ 2 T ||6||^ iterations. 
To complete the proof, we note that 



A{s°r i i m \\ < 



s° ■ 

mm 



□ 



C.4 Calls to the Solver 

In each call to Unshif t, we solve one system in a matrix of the form 

AS~ 2 A T = AS~ 2 A T + ms- 2 bb T 

— — LfLtp 

and in each call to Shift, we solve one system in a matrix of the form 

A§- 2 A T = AS~ 2 A T + ms- 2 bb T 

At the end of the interior-point algorithm we have one final call of the latter form. 

In order to say something about the condition number of the above matrices, we must bound 
the slack vector s. We are given an upper bound of T on the elements of s, so it remains to prove 
a lower bound: 

Lemma C.18. Throughout the InteriorPoint algorithm, s > 4gw ^ T ^ g 
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Proof. At all times duing the algorithm, we know from Lemma IC4I that the elements of s are 
bounded by a constant factor from the slacks at the current central path point s. In particular, 
taking into account Lemmas IC. 81 and IC. 141 we surely have s > ^s. So let us bound from below the 
elements of s. 

During the FindCentralPath subroutine, as we decrease z and expand the polytope ^l^ z , 
clearly the initial point s° remains in the interior of f2^ z throughout. Thus, by Lemma IC.31 we 
have s > ^s°, and so s > > -^s Q . 

Unfortunately, during the main part of the algorithm, as we increase z and shrink the polytope 
^6 2' ^ ne hiitial point may not remain inside the polytope. In particular, once we have z > b T y°, 
the initial point is no longer in £1® z , but we may define a related point (y z , s z , s z gap ) that is in ' z . 

Given our current point (y, s, s gap ) 6 ft^ z for z > b T y°, let us define r = -—j y ~% and note 

that < r < ^. We then define 



y z = ry u + (1 - r)y 



' s z ' 




c - A T y z ~ 




l°gap_ 




b T y z — z 





rs + (1 — r)s 
l 

2 



Hb T V-z) 



> 



Therefore Lemma IC.3I gives 



We then find 



1 . rs + (1 — r)s r r 

s > s = — > s l 

2m 2m 2m 



?gap . °gap 



2(b T y - b T y°) ~ 4nTU ~ 24nTU 

The last inequality follows because, when s gap decreased below | on the final step, using Lemma 
IC.4I we find that it certainly could not have decreased by more than a factor of ^ . 

We conclude s > \s > ^s° > ^±^8° □ 

We may now summarize the calls to the solver as follows: 

Theorem C.19. The InteriorPoint(^4, 6, c, y°, e) algorithm makes O ( ^fm log TU ™ — ) ca ll s 
to the approximate solver, of the form 

Solve (AS~ 2 A T + vv T , •, e(m~ 1/2 

and one call of the form 



Solve I AS'^A 1 + vv 1 , •, ft 



-2 aT , „..T r. I S min e 



m l / 2 n 2 T 2 U 2 , 

where S is a positive diagonal matrix with condition number O (^ T2u J n ' ) ; a nd v satisfies 

'u(mn) 1 / 2 \ 



O 



s° ■ e 

mm 



Proof. From Lemmas lC.17l and lC.12l the total number of solves is O ( y/m ( log TU ™ 1- log 



where we know from the FindCentralPath algorithm that s^ ap = 40 T ™J^ = O ( TU ™j£ /2 



mzn 
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As we noted above, all solves are in matrices that take the form AS 2 A T + vv T , where 
v = m^s-^b or v = m 1 / 2 s^S )- 1 ^ 

We know that s gap = 0(e) and s gap = Q(m), so we obtain the respective bounds 
_ c ( U(mn)V 2 \ / U{mn) l ' 2 \ 

\ e J \ S min J 

The condition number of S comes from Lemma IC. 181 and the upper bound of T on the slacks. 
The error parameter for the solver is (m^ 1 / 2 ) from the all NewtonStep calls. In the final 

solve, the error parameter is Smi ^ n — > yjj^ • 4SnmTU ' a § am invoking Lemma l"C. 181 □ 
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